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A  Fourier  representation  of  the  perturbation  effects 
of  the  moon's  eccentricity  on  a  satellite  in  orbit  about  L? 
was  constructed  and  verified  numerically.  It  was  then  incor¬ 
porated  in  a  modal  control  scheme  for  linear  systems  that  are 
time  periodic  which  was  developed  by  Shelton  (Ref.  f).  The 
feasibility  of  this  modified  control  scheme  was  then  verified 
by  computer  simulation  by  controlling  the  actual  non-linear 
motions  of  the  satellite  with  the  modified  controller. 
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Background 

Space  is  fast  becoming  a  part  of  our  everyday  lives. 
Sensor-bearina  satellites  placed  in  orbit  around  the  earth 
play  a  very  important  role  in  the  world  we  live.  From  ter¬ 
rain  and  resource  mapping  to  communication  and  strategic 
r econna i sance  our  lives  are  all  affected  by  our  evro’oi^s  in 
space . 

A  desire  to  use  the  minimum  number  of  satellites  to 
satisfactorily  accompl  ish  mission  objectives  and  maintain  the 
orbits  of  these  satellites  has  led  to  much  research  attemp¬ 
ting  to  develop  better  ways  to  accomplish  this  task.  A  num¬ 
ber  of  space  operations  require  total  coveraae  of  the  earth's 
surface  by  a  particular  set  of  satellites.  Currently,  sucb 
space  operations  are  conducted  by  satellites  positioned  in 
earth-synchronous  orbits,  but  not  only  does  this  not  provide 
total  global  coverage,  these  orbits  are  unstable  and  require 
some  sort  oF  orbit  maintaining  device  to  be  carried  by  the 
satellites.  This  control  system,  and  the  fuel  to  power  it, 
must  be  carried  as  an  integral  part  of  the  satellite,  taking 
away ,  pound  for  pound,  available  payload  weight  allocated  to 
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the  sensors  and  support  equipment  required  to  perform  the 
satellite's  primary  mission.  Also,  when  the  control  system's 
fuel  supply  is  exhausted,  the  sate1 lites  service  life  is  at 
an  end.  These  factors  combine  to  force  lesion  tradeoffs 
between  service  life  versus  mission  capability. 

In  recent  years  global  strategic  reconna i sance  and 
communication  have  become  increasingly  important  an!,  at  the 
same  time,  more  vulnerable  to  disruption.  With  the  Soviet 
Union's  increased  research  and  development  efforts  aimed  at 
developing  "sate1lite  killers"  a  less  vulnerable  position  for 
these  satellites  is  desirable. 

The  search  for  new  orbits  which  would  provide  a  max¬ 
imum  amount  of,  or  total,  earth  coverage,  hold  down  the  con¬ 
trol  cost  of  maintaining  these  orbits  and  still  decrease  vu1- 
nerability  to  attack  has  led  to  extensive  research  into  the 
"Lagrange"  points.  The  five  Lagrange  points  of  equilibrium 
shown  in  Figure  1  represent  theoretical  points  in  space,  in  a 
coordinate  system  that  rotates  with  the  earth-moon  system, 
where  the  combined  gravitational  effects  of  the  earth  and 
moon  will  keep  a  satellite  placed  at  one  of  these  points  at 
that  point.  However,  only  L^  and  L5  are  stable  and  wou’d 
provide  long  term  stable  orbits.  Since  the  Lanrange  points 
lie  in  the  earth-moon  plane  and  not  in  the  earth's  equitorial 
plane,  one  satellite  placed  at  each  of  the  points  L1,  Ld,  and 
L5  would  provide  "global"  coverage,  including  the 
strategically  valuable  North  Pole,  at  all  times. 
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Since  the  analysis  from  which  the  Lagrangian  points 
are  derived  is  based  on  the  assumptions  of  circular,  planar 
orbits  and  does  not  include  the  sun,  when  the  "restricted 
four  body  problem"  is  analyzed  all  of  the  Lagrangian  points 
become  unstable.  However,  the  points  do  provide  a  place  to 
start . 


Problem  and  Scope 

Shelton  developed  a  modal  control  scheme  for  linear 
systems  that  are  time  periodic.  This  scheme  was  then  applied 
to  a  satellite  in  orbit  about  the  earth-moon  Lagrange  point, 

L?.  However,  his  results  are  only  valid  for  a  set  of 

i 

dynamics  based  upon  his  specific  periodic  orbit.  The  free  j 

3 

lunar  eccentricity,  the  lunar  inclination,  and  the  eccen-  ! 

) 

tricity  of  the  sun  were  not  included  in  his  analysis. 

Since  the  free  lunar  eccentricity  is  the  major  con¬ 
tributor  of  in-plane  perturbations  about  the  periodic  orbit, 
this  study  will  focus  on  its  effects.  with  the  active  con¬ 
troller  developed  by  Shelton  operating,  the  periodic  orbit 
becomes  a  stable  reference  solution  and  classical  perturba¬ 
tion  techniques  can  be  applied  to  the  addition  of  the  free 
lunar  eccentricity.  If  the  solution,  including  the  first 
order  perturbation  terms,  is  well  behaved,  free  of  resonances 
and  unstable  roots,  then  the  controller  will  be  modified  to 
suppress  only  the  "free"  components  causing  the  original  in¬ 
stability.  Higher  order  terms  will  hopefully  be  several 
orders  of  magnitude  smaller  than  the  first  order  terms 
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leavinq  us  to  believe  if  the  controller  operates  satisfac¬ 
torily  in  the  lowest  order  perturbed  system  that  the  higher 
order  terms  are  negligible.  Applying  this  to  Shelton's 
system  we  will  again  arrive  at  a  control  system  with  effec¬ 
tively  zero  long-term  stationkeeping  costs,  but  that  is  valid 
in  the  "real"  representation  of  the  dynamics  and  not  just  in 
the  idealized  system. 


1 1 .  prohl. em  Anal  ys  i  3 

Model ino  Assumptions 

Shelton's  research  was  based  on  a  model  of  the  three 
attracting  bodies  (earth,  moon,  and  suni  that  placed  the  moon 
in  an  orbit  about  the  earth  that  is  periodic  but  not  cir¬ 
cular.  The  earth  orbit  around  the  sun  was  modeled  as  cir¬ 
cular  and  other  inclination  and  eccentricity  effects  were 
initially  ignored,  restricting  the  motions  of  the  four  bodies 
to  a  plane.  This  study  used  these  assumptions  to  rederive 
the  Shelton  controller  followed  by  a  perturbation  analysis. 

Coordinate  System 

Due  to  the  nature  of  the  periodicity  of  the  orbit 
about  L? ,  the  coordinate  system  used  in  Shelton's  research, 
and  in  this  study,  rotates  with  the  average  angular  velocity 
of  the  moon  and  the  center  of  the  system  is  located  at  the 
earth-moon  center  of  mass.  This  frame  allows  easy  visualiza¬ 
tion  of  the  periodic  nature  of  the  orbit  about  LI  (Ref.  f )  . 

Derivation  o f  the  Hami l tonian  (Ref.  G) 

In  deriving  the  dynamics  for  the  restricted  four-body 
problem,  the  three  attracting  bodies  were  assumed  to  have 
their  actual  masses  while  the  satellite  was  assumed  to  be 
massless.  This  assures  that  the  attracting  bodies  affect  the 
motion  of  the  satellite,  but  the  satellite  has  no  effect  on 
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the  motion  of  the  attracting  bodies.  The  geometry  of  this 
system  is  depicted  in  Figure  7. 

In  order  to  derive  the  Hamiltonian  the  kinetic  and 
potential  energies  must  first  be  constructed.  The  inertial 
point  in  this  derivation  is  assumed  to  be  the  earth-moon-sun 
center  of  mass.  This  results  in  the  inertial  position  vector 
of  the  satellite  to  be 

m. 


'7  _ 

rSat  =  R  -  -  P 


TT 


<D 


where  M  =  m^  +  (sum  of  the  masses  of  the  earth,  moon, 
and  sun,  respectively)  and  the  factor  m^/M  which  locates  the 
inertial  point  of  the  system.  Taking  the  inertial  time 
derivative  of  the  position  vector  in  the  rotating  frame 
yields 


ij  _  _  r^  _  _  m7r_i  _ 

dt,rSat>=''Sat=R  +(“  *  r>  +  ,T>  5  +  “  *  1-jH  T  I2) 

nt it • 

where  -—  "0  +  aj  x  (-— )  "p  is  the  transport  velocity  of  the 

M  M 

earth-moon  center  of  mass.  The  symbol  "Hi  denotes  the  angular 
velocity  vector  of  the  rotating  frame,  which  is  the  moon's 
mean  inertial  angular  velocity  vector,  and  the  superscript  r 
represents  differentiation  in  the  rotating  frame.  The 
kinetic  energy,  T,  of  the  satellite  is  defined  as 

T  =  1/2  m  (VSat  •  VS8t)  <’) 


P 


Because  the  satellite  is  assumed  to  be  massless,  the 


IT?  -  o 


IR  -  ^1  r! 
m  ^  +m  ^ 


IR  +  m2  r  1 

m  j  +m  2 


as  the  potential  energy  is  not  dependent  on  the  time  deriva¬ 
tives  of  the  coordinates,  the  momenta,  p^ ,  of  the  Hamiltonian 
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where  q.  is  the  time  derivative  of  the  i  coordinate.  The 
coordinates  used  were  the  x,  y,  and  z  components  of  the 
satellite's  position  as  depicted  in  Figure  2.  In  this  frame 
the  momenta,  p. ,  are  the  actual  inertial  velocity  components, 
a  fact  that  Shelton  used  to  derive  his  controller. 

The  Hamiltonian  is  given  by 


H  =  Z  Piqi  -  Tg  +  Vc 
i  =  1 


Combining  equations  M),  (*)  and  (7)  gives  us  the  Hamiltonian. 


^  9  T  . 

H  =  ~(px  +py  +pz  }  +  “  fpxRy  "  pyV  +^r  (  Px"*0  y)px 


+  f  Py+^jfJPy1 


IT?  -  P  I  IR  -  ml  r  I  I R  +  n2  r  I 


m  j  +m  j 


m  2+m  p 


From  this  Hamiltonian  the  equations  of  motion  of  the 


satellite  were  formed. 

Equa t ions  of  Motion 

The  equations  of  motion  for  the  satellite  are  derived 
from  the  Hamiltonian  by  forming 

.  3H 

qm=  —  m=l ,n  (10) 


.  3  H 

pm  =  _  m-lf? . n  (11) 

3qm 

where  n  represents  the  number  of  generalized  coordinates 
describing  the  system  (Ref.  5:94).  These  equations  become 

x  (t)  =  'f  (x(t)  ,t)  (12) 

when  put  in  state  variable  form.  The  state  vector  x 
is  made  up  of  the  coordinates  and  momenta,  q's  and  p's, 
and  7  contains  the  partial  derivatives  of  the  Hamiltonian 
defined  by  (10)  and  (11).  In  this  restricted  four-body 
problem  the  Hamiltonian  is  an  explicit  function  of  time  so 
the  equations  of  motion  are  explicit  functions  of  time  also. 

Shelton  assumed  the  existence  of  a  periodic  solution 
to  (1?),  Xp(t).  Defining 

6x (t)  =  x (t)  -  xn (t)  (13) 
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and  differentiating  (1?)  to  get 


Sx(t)  =  x(t)  -  xp(t)  (1^1 

Substituting  f  1 3 >  and  (14)  into  (12)  results  in 


•  • 

x0(t)  +  Sx(t)  *  f(x0(t)  +  6x(t),t)  (15) 

Expanding  this  equation  in  a  Taylor  Series  about  j?p(t) 
produces 

1.  —  —  —  3  ^  i  _ 

X0(t)  +  <5X(t)  =  f  ( X  p  ( t )  ft)  +  |<$  X  ( t )  +  H.O.T.  (15) 

3  x  ( t )  xp(t) 


but  from  (12)  we  see  that 


x0  (t)  =  T (xp  (t)  ,t) 
which  reduces  (15)  to 


(37) 


5  X  (t)  a  Aft) 5  x  (t)  (IP) 

including  only  first  order  terms,  with 

3f  | 

Aft)  =  —  I  (19) 

3X(t)  x0(t) 


Displacing  the  periodic  solution  ,  xp(t),  at  tp  an 
amount,  6x,  and  expanding  the  solution  in  a  Taylor  series 

about  x ^  C t )  yields 

„  3xft)  . 

X(t)  =  Xpftl  +  -  I  5x(tp)  +  H.O.T.  (20) 

3x  ftp)  xp(t) 

Truncating  to  first  order  and  using  n 31  produces 

«x(t)  =  *  (t,tp)  6x(t0)  (21) 

where 


3x(t)  , 

Mt,tn)  =  -  j  (22) 

3x(tp)  xp<t) 

The  matrix,  $,  is  called  the  state  transition  matrix  and  is  a 
function  of  only  the  initial  and  final  times  being  con¬ 
sidered.  The  properties  of  the  <i>  matrix  are 

$(t,tp^  =  ♦ft^.j)  «{tj,tpJ  (27) 

and 

*(t0,tp)  =  i  (?n 


where  I  is  the  Identity  matrix 
respect  to  time  produces 


Differentiating  (21)  with 


5xft)  =  5  x  ftp)  {??) 

as  Sxft^)  is  a  constant. 

Floquet  Theory 

The  four-body  problem  Hamiltonian,  derived  in  a 
previous  section,  is  an  explicit  function  of  time  and  not  a 
constant  of  the  motion.  Enerqy  is  not  conserved  and  no  other 
constants  of  the  motion  exist  fRef.  5:41$).  The  equations 

of  motion  that  are  derived  from  the  Hamiltonian  are  also  time 
varying.  The  variational  equations  of  motion  are  written  in 
state  vector  form,  derived  previously,  as 

6xft)  =  Aft)  6xft)  (?$) 

This  system  is  periodic  and  the  periodicity  is  contained  in 
the  coefficient  matrix  A(t)  where  t  is  defined  as  the  period 
of  the  system.  Differential  equations  such  as  (?$),  that 
have  periodic  coefficients,  can  be  reduced  to  the  case  of 
constant  coefficient,  linear  differential  equations,  by  the 
theory  of  Floquet  (Ref.  ?:^0).  Floquet's  theory  states  the 
general  solution  to  (?$)  is  of  the  form 

n  _  At 

6x  ft  )  =  l  8  .E.  ft )  e  3  {?’>) 

j  =  l  D  3 

where  8^  and  A.,  are  constants  and  the  E^  are  periodic  func- 
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tions  with  the  sane  period,  t.  The  n  can  he  arranger"  in  a 
matrix  of  the  Jordan  canonical  form  which  basically 
diagonalizes  the  matrix  (Ref.  f:?^'7).  The  constant  A ...  is 
called  a  character istic  or  Poincare  exponent  of  Aft). 


needed 

Letting 


In  order  to  compute  E_.  and  Shelton  assumed  that  he 
to  excite  one  mode  of  the  system  described  by  (?f). 
8^=0  for  j^i  and  8  .=1  yields,  from 


5  x  ( t ) 


Ej  (t)e 


Ait 


(  2?) 


From  (21),  with  tp=0 


fix (t)  =  i(t ,0)Sx(0)  (??) 


and  equating  (2P)  and  (29),  yields 


-  \t 

E.  ( t )  e  =  f(t,0)fi  x  ( 0 )  (3d) 


-  $(t,0)  E.  (0) 


(31) 


Since  E^(t)  is  a  periodic  vector,  E^(x) 

A  .T 

5  x  ( x  )  =  E .  (  t)  e  1 

l 

X  .x 

=  E.  (0)  e  1  =  «t,0)E.  (0) 


E.(d),  results  in 


(3  2) 

(33) 


Rearranging  (33)  yields 


r  $(T  ,  0)  -fe 


0 


•  »  ' 


f 


xiT 

nTFjio) 


<  3*) 


which  is  an  eigenvalue  problem  with  the  eigenvalues  of  4>(t,P) 
being  e  .  The  eigenvectors  of  $(xr^)  are,  therefore,  the 


vectors  E^ (0)  . 

To  complete  the  solution,  E^ft)  must  be  computed  over 

one  period.  Since  "E.  ft)  is  periodic  this  determines  E.  ft) 
for  all  time. 

Combining  equations  (?5)  and  (?.P)  produces 


d  _  A  { t  _  t 

—  (E.  me  )  =  Aft)  ( E .  ( t )  e  )  (35) 

dt 


Defining  A( t )  as 


Aft)  =  f£1 ft)  I E? ft)  I  .  .  . ! En  ft)  1 

and  incorporating  Aft)  into  (95)  yields 

.  Jt  Jt  Jt 

A(t)e  +  Aft)Je  =  A  f  t )  A  (tie 

and  reducing  further  gives  the  d ; f f e ren t  i  a  1  equation 


(3*) 


f  3"7) 


A  ft)  +  Aft) J  =  Aft) A  ft)  (3P) 

Shelton  used  the  Floquet  theory  described  above  to 
investigate  the  stability  of  the  orbit.  His  work  demon¬ 
strated,  as  did  work  previously  done  by  Wiesel  and  Smith, 
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that  the  orbit  was  unstable.  One  of  the  system's  Poincare 
exponents  is  a  positive  rea]  root,  while  its  conjugate  is 


negative  and  the  others  are  purely  imaginary. 

Table  I.  Poincare  Exponents  of  Shelton's  Orbit 


Mode 

Poincare  Exponen 

Planar 

0.0  -  l.C392i 

0.0  +  1 . 03??i 

Planar 

2. 3921  +  0.  Oi 

-2.3921  +  O.Oi 

Out  of 

0.0  -  1.12471 

Plane 

0.0+  1.12471 

In  Shelton's  orbit,  the  unstable  exponent  was  designated  X, 
and  is  the  root  that  his  controller  affected  (Ref.  *) .  This 
gave  the  following  equation  when  the  controller  was  incor¬ 
porated 


«x(t)  =  A(t)«x(t)  +  B(t)uft)  (39) 

where  u(t)  is  the  control  vector  and  B(t)  is  used  to  desig¬ 
nate  to  which  states  the  control  will  be  applied.  In  order 
to  design  the  control  system  so  only  the  unstable  mode  is 
affected  by  the  controller,  a  new  variable  was  defined. 

-1 

net)  =  A  (t)  x(t)  (*0) 
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where  nft)  is  the  modal  vector  and  Aft)  is  as  defined  in 
(36).  Taking  the  time  derivative  of  equation  MO)  yields 

•  •  • 

S  x  ( t )  =  \  (t)n(t)  +  AftJnft)  (41) 

Equating  equations  (39)  and  Ml)  and  solving  for  "nft),  to  get 
a  differential  equation  in  terms  of  the  modal  variables, 
results  in 

^  -1  .  -1 

nft)  =  A  (t)  (A(t)A  (t)  -A(t)lnft)  +A  (t)Bft)u(t)  M?) 
From  (3?) 

* 

Aft)  =  Aft)  Aft)  -  Aft)  J  M3) 

substituting  M?)  in  M?)  produces 

^  _  -1 

n  ft)  =  J n ( t )  +  A  (t)B (t)u  ft) 

where  J  is  the  constant  diagonal  matrix  of  Poincare  ex¬ 
ponents.  Since  the  system  is  diagonalized,  feeding  back  a 
sinqle  unstable  mode  can  be  accomplished,  making  it 
relatively  simple  to  implement  the  control. 

Shelton  implemented  this  control  scheme  by  selecting 
uft)  =  U  =  kn^  where  k  is  the  constant  control  gain.  This 
gave  the  following  differential  equation 


18 


controlled  system,  a  perturbation  analysis  was  performed  on 

the  system  and  the  controller  was  modified  to  ignore  the 
oscillations  forced  on  the  orbit  by  the  moon's  eccentricity. 
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i 


III. 


Perturbation  Analysis 


The  solution  of  a  differential  equation  of  the  sort 


x  =  A  x  +  F  (*5) 

contains  a  free  (homogeneous)  and  a  forced  (particular)  solu¬ 
tion  such  that  the  total  solution  is  of  the  form 


X  =  X.  +  x_ 
h  p 


(46) 


with  Xp  being  the  particular  solution  associated  with  a 
particular  F  (forcing  function).  Shelton's  controller  is 
designed  to  suppress  both  the  free  and  forced  components  of 
the  unstable  mode.  The  controller  was  designed  for  the 
idealized  case  with  the  moon  in  a  circular  orbit.  when 
Shelton's  controller  is  implemented  with  the  moon  in  its 
real,  slightly  eccentric,  orbit  the  controller  does  not  keep 
the  satellite's  orbit  stable.  In  order  to  incorporate  the 
perturbations  caused  by  the  moon,  in  the  variational  equa¬ 
tions  the  lowest  order  perturbation  terms  evaluated  on  the 
periodic  orbit  will  be  included  such  that 


.t  -i 

5x  =  fflft)  +  BK  A  (t)H  x  +  P(t)  (-37) 

where  K  =  f  0  0  k  01  when  only  0,  is  fed  back, 

P(t)  represents  the  perturbing  effects  of  the  moon's  eccen- 


2 


fci 


tricity  and  K A  ( t )  x  is  the  Shelton  controller.  Defining 


A ’  (t)  =  Aft)  +BKA  ft) 


and  applying  the  Floquet  transform  to  the  modal  variables  for 
this  controlled  system  (see  section  on  Floquet  theory),  we 
find 


n'  It)  =  j'n 1  (t)  +  A'  ft)  Pft) 


where  the  primes  refer  to  the  closed  loop  system,  for  a  given 
gain  k.  This  reduces  the  system  (^7)  to  constant  coefficient 
form  (Ref .  9) . 

The  solution  to  the  constant  coefficient,  linear  dif¬ 
ferential  equation  is  of  the  form  f^l'.  The  forcing  function 


Pft)  was  developed,  to  lowest  order,  using 


Pft)  =xR(t)  -  xpp  ft) 


where  xR(t)  is  the  derivative  of  the  solution  to  the 
variational  equation  with  the  noon's  actual  eccentricity 

included  and  Xp^ft)  is  the  derivative  of  the  solution  with 
the  moon  in  its  idealized,  circular  orbit. 


A'  Pft)  was  obtained  by  evaluating  (50)  at  thirty 
equally  spaced  time  steps  around  the  orbit  while  holding  the 
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angle  associated  with  the  phase  of  the  moon  constant, 

-J 

multiplying  times  A'  (t) ,  and  then  constructing  a  Fourier 
series  to  represent  the  elements  at  each  time  step.  Keeping 
only  the  first  N  terms  of  the  series  yields 


-1  _  N-l  N-l_ 

A ' (t)  P(t)  =  l  T.'  .cos  (a,)  +  E  B ' . sin (a ■ )  (51) 

i =0  i =0 


where 

a.  =  iw0t 

This  procedure  was  repeated  until  the  Fourier  series  for  a 
number  of  small,  constant  step  increments  of  the  phase  of  the 
moon  from  0  to  radians  were  obtained.  Another  Fourier 

series  was  then  constructed  from  the  coefficients,  A^'s  and 
B  r ' s ,  of  these  Fourier  series.  The  result  was  a  double  sum 
Fourier  series  represented  by 

-1  _  N-l  M— 1 

A.ft)P(t)  -  £  5:nrAnmcosfan)cos^bm)  +  Bnmcos  sin  (V 
n=0  m=P 

+  C  sin(a„)cos(b  )  +  D  sin  (a„ )  sin  (b_)1  f  5  2 ) 

nm  n  m  n  mJ 

where  the  coefficients  A  ,  B  ,  C  ,  and  D  are  the  coef- 

nm  nm  nm'  nm 

ficients  resulting  from  the  construction  of  the  overall 

Fourier  series  (52).  The  terms  a„  and  b  are  associated  with 

n  m 

the  phase  of  the  orbit  and  the  phase  of  the  moon,  respec¬ 
tively,  and  are  defined  as 
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(53) 


f  * 


an  =  nw0t 

bm  "  m(wEfc  +  °10>  f54) 

where  n  and  m  are  the  summation  indices,  w^  is  the  periodic 
orbit  frequency,  w^  is  the  frequency  of  the  line  of  apsides 
of  the  lunar  orbit,  and  is  a  constant  phase. 

Since  the  product  of  two  trigonometric  functions  can 
also  be  represented  as  the  sum  and  difference  of  two  trig 
functions  the  following  relationships  were  used  to  expand 
-1. 

A ' ( t )  P(t). 


1  1 

sinfa'sin(b)  =  —  cos(a-b)-  —  cos  (a+b) 
2  2 


(55) 


1  1 

cosfa)cos(b)  =  —  cos(a-b)  +  —  cos(a+b) 
2  2 


(55) 


1  1 

sin(a)cosfb)  =  —  sin(a+b)  +  —  sin(a-b) 
2  2 


(57) 


1  1 

cos(a)sin(b)  =  —  sin(a+b)  +  —  sin(a-b) 
2  2 


(5P) 


Defining 


,nm  =  a  +  b  =  (nw0  +  mwE)t  +  m010 


(59) 


Sm  s  a  -  b  =  (nw0  -  "wE)fc  -  m010 

and  substituting  (55)  through  (50)  into  (52)  produces 


(50) 
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'♦"I"  1  .f 


-1 


N-]  M-l 


A*  (t)  P(t'=l/2  E  E  r  _'cos((i  )  +  (A  -C  )  cos  ( <(>_  ) 

-  .  nm  nm  run  nm  run  nm 

n  =  P  m=0 


+ (B  -D_m) sin  (  i>  )  +  (B +D  )  si  i  ( <t>  )  ) 

nm  nm  nm  nm  nm  nm 


(6] ) 


From  (49)  we  see  that 


-1 


n ’ (t)  =  J • n •  ft)  +  A» (t) P(t) 


( £2) 


Since  J'  is  of  the  form  (Ref.  5) 


J'  = 


0  a1  o  o 
a2  0  n  0 
o  o  cx^  n 
0  n  0  a, 


( ^3) 


where  the  a's  are  the  poincare  exponents  of  the  controlled 
system.  Since  n-,(t)  is  used  by  the  controller  solving  (^>2) 


for  n^p(t)  results  in 


R  3p  ( t )  =  «?n  +  fA«  (t)P(t) 1 3 


(<S4) 


Assuming  a  solution  of  the  same  type  as  the  forcing  function 
and  takina  its  derivative  we  obtain 


NM 


n-,p(t)  =  E  rwLcos  (<j>L) +VLsin  f  <|>L) +YLcos  (ij/L) +ZLsin  (ij>L 
L=  1 


)  1 


(«?) 
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.  NM 

n-,pft)  *  I  r-WL4)  Lsin  L) +VL(j>  Lcos  (<{>  L) 


«  • 

Lsin  L) +ZL^  Lcos  L)  ]  (33) 

where  WL,  VL,  YL,  and  ZL  are  the  coefficients  we  are  solving 
for  and  L  is  related  to  the  nm  designator  in  (31). 

Substituting  (31),  (35)  and  (33)  into  (3*)  and 

equating  like  terms  yields  (  L  corresponding  to  a  particular 
nm,  L  =  Mn  +  m  +  1) 


sin  (4>. )  : 

-WL»L 

=  a ,V.  +  1/2  (B  +D  ) 

?  L  '  nm  nm' 

(37) 

cos  (<p  L)  : 

VL»L  - 

a-,Wr  +  1/2  (A  -C  ) 

3  L  nm  nm 

( 3P ) 

sin  (^ L )  : 

-Vl 

=  a  zr  +  1/2 (B  -P  ) 

?.  L  nm  nm 

(39) 

cos (^L) : 

2L*L  * 

a,Yr  +  1/2 (A  +C  ) 

?  L  '  nm  nm 

(70) 

Solving  (37)  through  (70)  for  W^,  VL,  Y^  and  ZL  and  sub¬ 
stituting  in 


'P 


L 


L 


we  obtain 


=  nWp  +  mw£ 

=  nWp  -  mWg 

the  coefficients  for  the 


particular 


(71 ) 

(72) 

solution  of 


n3 


2  3 


(t)  as  follows 


(73) 


WL  " 


-a  V--1/2(B  +D  ) 
3  L  nm  nm 

(nw0  +  mwE) 


VL  - 


a.Wr+l/?(A-C) 
3  L  nm  nm 

(nwp  +  mv/g) 


(74) 


yl  - 


-a  Zr -1/7 ( B  -D  ) 
3  L  nm  nm 

(nwp  -  mwE) 


ZL  - 


a,Yr+l/2(A  +C  ) 
3  L  nm  nm 

(nw0  -  mwg) 


(75) 


(75) 


Numerically  evaluating  (73)  through  (75)  ?nd  using  them  with 
(55)  gives  us  a  solution  to  n3  (t).  To  insure  that  this 
solution  was  correct  n3p(t1)  was  numerically  calculated  and 
using  the  approximation 


^3p(tl] 


^  3p  ( 1 1 )  ~  a-jpftj+At) 

At 


(77) 


and  substituting  these  values  into 

P3(t1)  =  A  ( t x )  fn 3  ( t x )  -  a,n3p(t])l  (7?) 

This  calculated  value  was  compared  with  the  value  for  P7(t) 
calculated  from  (50)  .  The  numbers  matched  to  six  decimal 
places  using  a  At  =  lO-^  sec.  This  makes  me  resonably  sure 
that  the  correct  solution  for  n,  ft.)  was  formulated.  The 


derivations  of  n^pft)'  nppft)  anrf  n^p(t)  are  contained 
Appendix  A. 


IV.  Implementation  of  the  -Mod i  f  i  ed  Controller 


Having  solved  for  the  particular  solution  due  to  the 
perturbation  effects  of  the  moon  Shelton's  controller  was 
modified  to  incorporate  this  particular  solution.  Shelton's 
controller  stabilizes  the  idealized  orbit  by  zeroing  out  the 
difference  between  the  satellite's  true  orbit  and  the  stable 
controlled  periodic  orbit.  When  the  moon's  perturbation  ef¬ 
fects  are  added  in,  the  controlled  orbit  is  no  longer 
stabilized.  The  Fourier  analysis  of  the  perturbinq  force 
showed  it  to  be  purely  oscillatory,  with  no  secular  terms,  to 
first  order.  If  we  could  calculate  these  terms  and  subtract 
them  from  the  controlled  term  we  could  stabilize  the  orbit 
with  the  moon's  proper  dynamics  involved. 

The  Shelton  controller  stabilizes  the  idealized 
system  by  calculating  the  control  term  given  in  MO 

-1 

CT.  =  k  A  ft)  Bn.,  ft)  f 79) 

which,  with  the  moon  in  its  actual  orbit,  yield 

-1 


CT.  =  kA  (t)Bfn?Hft)  +  n,  (t)l 

fPO) 

driving 

the  system  unstable.  Having 

solved  for 

the 

par- 

t icul a  r 

solution  f65)  the  controller  was 

impl emented 

by 

sub- 

tracting  off  n^pft)  which  yields 

2R 


.  1  A . 


mm 


7 


-1 

CT.  =  kA  (t)Bn ?H(t)  (PI) 

This  system  produces  a  new  stabilized  orbit  which  is  not  the 
idealized  periodic  orbit  that  Shelton  derived. 
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the  controller  both  the  actual  integrated  solution  with  the 


modified  controller  operating  and  the  particular  solution, 

n  (t)  ,  were  plotted  to  see  if  £he  integrated  solution  would 
track  the  calculated  particular  solution  (Figure  3).  All 
four  modal  variables  were  plotted  to  insure  that  the  cal¬ 
culated  solution  was  indeed  tracking  the  integrated  solution 
for  each  variable-  As  can  be  seen  in  Figure  3  the  calculated 
and  integrated  solutions  match  very  well. 

During  the  sixth  orbit  the  integrated  solution 
becomes  unstable  (Figure  ^).  This  occurs  due  to  a  combina¬ 
tion  of  two  factors,  the  gain  selected  for  implementation  of 
Shelton's  controller,  k  =  .<*59,  was  small  and  the  particular 
solution  was  larger  than  expected.  Because  the  particular 
solution  was  so  much  larger  than  expected,  the  error  in  the 
calculated  term  is  also  quite  appreciable.  This  is  because 
the  perturbation  effects  were  assumed  to  be  first  order  and 
were  evaluated  with  respect  to  the  periodic  orbit,  which 
would  have  been  a  good  assumption  if  the  particular  solution 


Figure  3(a)  ETA  vs  Time 


were  two  orders  of  maqnitude  smaller.  A  phase  portrait  for 
the  Shelton  Controller  with  a  gain  of  0.5  is  depicted  in 
Figure  5  (Ref  5).  As  can  be  seen  in  this  figure  some  of  the 
trajectories,  on  the  order  of  =n  ^,  =  1 0-^ ,  do  not  converge  to 
the  origin.  The  combination  of  factors  mentioned  above  put 
this  solution  on  one  of  the  trajectories  that  is  unstable. 

Due  to  the  combination  of  factors  described 
previously  the  system  is  not  stable  beyond  six  orbits.  This 
situation  could  possibly  have  been  avoided  if  the  particular 
solution,  np3 (t) ,  used  in  the  controller  was  more  accurate  or 
if  a  larger  gain  was  selected  for  the  original  controller  at 
the  start  of  the  analysis.  Selecting  a  larger  gain  would 
have  increased  the  stable  operating  region  of  the  controller. 
However,  there  was  no  way  of  actually  knowing  ahead  of  time 
how  large  or  accurate  the  particular  solution  would  be. 

In  spite  of  the  failure  of  the  controller  to  operate 
beyond  six  orbits,  it  did  exhibit  excellent  control  charac¬ 
teristics  up  to  the  time  it  became  unstable.  In  Figure  5  we 
can  see  that  the  integrated  and  calculated  solutions  track 
very  well.  In  this  case  the  integrated  solution  was  started 
on  the  calculated  solution.  To  see  if  the  controller  would 
drive  the  integrated  solutior  to  the  calculated  solution,  I 
started  the  integrated  solution  off  of  the  calculated  solu¬ 
tion  and  the  controller  did  drive  it  back  to  the  calculated 
solution  (Figure  7).  Figure  P  shows  the  difference  between 


the  calculated  and  integrated  solution  for  five  orbits  and 


gure  6  ETA3  Started  on  Orbit 
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Figure  8(a)  Difference  Between  Calculated  and  Integrated  Solutions 


Figure  9  shows  the  amount  of  control  being  expended  to  main¬ 


tain  the  orbit.  The  integrated  control  gains  (ICG)  shown  in 
Figure  °  are  defined  as 

1  t 

ICG(t)  =  —  /Icontrol  acc.l  dy 
t  0 

which  is  equal  to  the  average  control  acceleration  required 
for  system  operation.  Expressing  these  values  in  terms  of 

_7 

g's,  yields  an  average  control  cost  of  approximately  10  g 
which  is  comparable  to  a  typical  earth-synchronous 

_7 

satellite's  control  cost  of  10  g  (Ref.  2:^3). 

Conclusions 

As  shown  in  Figures  7,  p  and  9  the  modfied  controller 
does  work  very  well  over  this  five  orbit  period.  Because  the 
initial  selection  of  gain  for  the  Shelton  controller  was  too 
low  I  did  not  achieve  the  long  term,  low  cost  solution  I  was 
expecting  to  develop.  However,  because  of  the  initial  well 
behaved  aspect  of  the  controller,  I  believe  the  method  is 
sound.  In  order  to  fully  develop  this  control  system  and 
prove  that  it  is  feasible,  I  recommend  that  a  follow  up  study 
be  done  using  a  higher  gain  and  also  include  the  effects  of 
the  inclination  of  the  moon’s  orbit. 


/») 


Figure  9(a)  Controlled  System  (Stalled  on 


f 
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Figure  9(b)  Controlled  System  (Started  off  Orbit) 
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Append i x  A 

Derivation  of  the  Particular  Solution 

In  deriving  the  particular  solution  to  (62)  we  can 
see  that  the  form  of  the  solutions  n^ft)  and  n^ft)  are  the 
same  except  for  a  change  of  indices.  If  I  change  the  3's  to 
4's  in  (64)  on  page  24,  and  follow  through  the  rest  of  the 
derivation  on  pages  24  through  26,  the  solutions  are 
generically  the  same. 

The  differential  equations  involving  n^-ft)  and 
n2pft)  are  coupled  and  of  the  form 

.  -1 

nlp(t)  =  “ln2p(t)  +  rA  •  (t)  Pft>  1  2  (P2) 

-1 

n2p(t)  =  «2T1lp(t)  +  rA’(t)P(t)1?  (P3) 

Taking  the  time  derivative  of  f P 2 )  and  substituting 

(P3)  produces 

n  ipft)  =  ala2nlp(t)  +  “lF2  +  F1  (p/>) 

-1 

where  is  the  time  derivative  of  rA'(t)P(t)lj  and 
-1 

F?  =  r  A ' f t ) P  ft  1  1 2  («5) 


Assuming  a  solution  to  (P*3)  of  the  form 


NM 

=  E  fWLcos  (<pL)  +VLsin  (<fiL)+YLcos  (ipL)  +ZLsin  fipL)  1  (8 6) 


taking  time  derivatives  of  (P6)  yields 


NM 

nlp(t)=  E  r-WL$Lsin  f  4.  L)  +VL(tLcos  f  4>  L1  -YL^Lsin  ( ^  L)  +ZL\|)Lcos  f  ^  L )  1  (97) 


. .  NM  .  2  .2  .2  .2 

nlp  (t)«  E  f-WL<))Lcos  (<})L)-VL$Lsin  (*L)  -VL*Lcos  (*L)  -Z^sin  (H»L)  ]  fPP) 


wi  th 

4>  =  (nw^  +  mwg)t  +  mPlc  (p9) 

1J1  *  (nwp  +  mw£)t  -  m01(3  (90) 

<j>  «  nwQ  +  mw£  (91) 

^  =  nwQ  -  mw£  (92) 

<j>  =  'P  =  0  (93) 


-1 

taking  the  time  derivative  of  T A' ( t) P (t ) 1 ^  ,  from  (61) 
page  2^,  yields 

'l  -  -  "  ;-*nm'An»l+cnml>sln','nm>-»nn'Anr,l-cnml>SInl*nm' 

2  n=0  m=0 


+  lb  (B  ,-D  ,)COS(<j>„  )  +  $ _ (B  ,  +D__  1  )  cos  (  ♦  )  1 

vnm  nml  nml'  nnr  nm  nml  nml  nm 


-1 

with  r  A' (t) P(t) 1 2  equal  to 

1  N-l  M-l 

F2  =  “  Z  1  !!fAnm2+Cnm2)cosf'J'nm)  +  fAnm2'Cnm2)cos((*>nm) 

2  n=0  m=0 

+  'Bnm2-Dnm2>si"'*nm'  +  ,B™?+Dnm2>SI"'V>'  «951 

taking  (87)  through  (95),  substituting  then  into  (8^)  and 
equating  like  terms  produces 


sin  ( <f>r  )  : 

,  1 

-VL(nw0+mwE)  =  «1«,VL+-c,1(Bnm2+Dnn|2)-rnwf,WE)  (Antll-Cnml)  («<>) 


COS  (  <j>L) 


2_ 


-WL(nwp+nwE)  =  a1a2WL+-a1(Anm2-Cnm2)  +  (nw0+niwEMBniTl]+Dnml)  (97) 


sin  (i(»L) 


2. 


_ZL(nwO+mwE'  ala2ZL+^al fBnm2_Dnm2) -(nwp~mwE) (Anmi+Cnm) 5  r?R) 


cos (^L) 


7  _ 


1 


-YL(nwp-mwE)  =  «  ?YL+-a  i  ^  Anm2+Cnm2^  +  f  (Bnml-Dnml  ^  (?9) 


solving  for  the  coefficients  of  the  particular  solution  from 
(95)  through  (99)  we  obtain 


W 


1L 


1/2  a 


l(Anm7-Cnm2>+  l/^fnwp+mwE) (Bnml+Dnml) 

- 

fa  2  +  (nwp+mwE)  1 


(100) 


V1L  " 

1/?  “ , (Bnm,+Dnn7)-  l/^(nw0+mwE) (Anml-Cnml) 

-  — - - - - ~ 

fa  ^a  2  +  (nw0+mwE)  1 

(101) 

1/2  a,(Anffl,+Cnffl,)+  l/?(nw0-mwE) (Bnnl-Dnnl) 

(102) 

II 

«— ( 

>< 

-  - - —  ~T.  ” 

raia2  +  (nw0-mwE)  1 

1/2  a,  (Bnm9-Dnm2)-  l/?(nwp+mwE)  (Anml+Cninl) 

(103) 

Z1L  * 

- - - -  -?~ 

fa  iCt  p  +  (nWp-mwE )  > 

Numerically  evaluating  (100)  through  (103)  and  using  them 

with  (PS)  gives  us  a  solution  to  nlp(t}.  Equations  (100) 

through  (103)  can  be  used  to  provide  the  coefficients  for 

n  (t)  by  switching  the  sub-indexes,  1  and  2. 

^P 


ip 
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